rm(list=ls())
library(bucky)
library(readstata13)
library(dplyr)
library(tidyr)
library(ggplot2)
library(data.table)
library(xtable)
library(marginaleffects)
location <- "/Users/dillonlaaker/Dropbox/shocks_replication/"
setwd(paste(location, "data/final", sep = ""))

load("ess_main.Rdata")

baseline <- " + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround)"


###################################
#############TABLE A1##############
###################################
##Generates table with description of independent (control) variables

Variables <- c("Female", "Minority", "Parental Citizenship", "Parental Education", "Secondary Educ.", "Avg. Education", "Educational Quality", "Ideology", "Income", "Education Years")
Descriptions <- c("Gender of respondent", "Minority status of the respondent", "Whether the respondent's parents are citizens of country","Sum of father's and mother's educations", "Whether the respondent completed at least 12 years of schooling", "Average years of education in the total population aged 15 years and older", "The extent that high quality basic education is guaranteed to all.", "Ideology of respondent", "Measures the respondent’s net total income", "Respondent's years of schooling")
Coding <- c("0-Male; 1-Female", "0-Not Minority; 1-Minority", "0-Parents Non-citizens; 1-Both Parents Citizens", "Range: 0(min)-8(max)", "0-Respondent completed less than 12 years of schooling; 1-Respondent completed at least 12 years of schooling", "Continuous Measure", "0 (high equality) to 5 (low equality)", "0-Left; 10-Right", "1-lowest income bracket; 10-highest income bracket", "Continuous variable: 0-56")

data <- data.frame(Variables, Descriptions, Coding)
table <- xtable(data, type = "latex", latex.environments = "left", caption = "")
x <- nrow(table) - 1
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE, hline.after=1:x, file=paste(location, "paper/tables/tableA1.tex", sep = ""))


###################################
#############TABLE A2##############
###################################  
#Generates table with summary statistics of variables

ess <- ess  %>%
  mutate(fhighskilled = ifelse(foccupation14=="High-skilled", 1, ifelse(!is.na(foccupation14), 0, NA))) %>%
  mutate(flowskilled = ifelse(foccupation14=="Low-skilled", 1, ifelse(!is.na(foccupation14), 0, NA))) %>%
  mutate(fmodskilled = ifelse(foccupation14=="Medium-skilled", 1, ifelse(!is.na(foccupation14), 0, NA))) %>%
  mutate(mhighskilled = ifelse(moccupation14=="High-skilled", 1, ifelse(!is.na(moccupation14), 0, NA))) %>%
  mutate(mlowskilled = ifelse(moccupation14=="Low-skilled", 1, ifelse(!is.na(moccupation14), 0, NA))) %>%
  mutate(mmodskilled = ifelse(moccupation14=="Medium-skilled", 1, ifelse(!is.na(moccupation14), 0, NA))) %>%
  as.data.frame()

ess <- select(ess, immigration, shock4_1825, shock4_1825comb, fhighskilled, fmodskilled, flowskilled, mhighskilled, mmodskilled, mlowskilled, female, minority, pcitizen, peduc, age, cohort, secondary, eduq_1017, aveedu_1017, educationy, income, lrideology, datalim)

ess$age[ess$age<26] <- NA
ess$datalim[ess$datalim==1] <- NA
ess <- ess[!is.na(ess$datalim),]
ess <- ess[!is.na(ess$age),]
ess$datalim <- NULL

nonNAs <- function(x) {
    as.vector(apply(x, 2, function(x) mean(x, na.rm=TRUE)))}
mean <- nonNAs(ess)
nonNAs <- function(x) {
    as.vector(apply(x, 2, function(x) sd(x, na.rm=TRUE)))}
sd <- nonNAs(ess)
nonNAs <- function(x) {
    as.vector(apply(x, 2, function(x) min(x, na.rm=TRUE)))}
min <- nonNAs(ess)
nonNAs <- function(x) {
    as.vector(apply(x, 2, function(x) max(x, na.rm=TRUE)))}
max <- nonNAs(ess)
nonNAs <- function(x) {
    as.vector(apply(x, 2, function(x) length(which(!is.na(x)))))}
length <- nonNAs(ess)
names <- c("Anti-Immigration Attitude", "Economic Shock (18-25)", "Economic Shock Count (18-25)", "Father High-Skilled", "Father Mod-Skilled", "Father Low-Skilled", "Mother High-Skilled", "Mother Mod-Skilled", "Mother Low-Skilled", "Female", "Minority", "Parental Citzenship", "Parental Education", "Age", "Cohort", "Education 12 years or more", "Educational Equality (10-17)", "Education Level (10-17)", "Education Years", "Income", "Ideology")
table <- data.frame(names, mean, sd, min, max, length)
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA2a.tex", sep = ""))

load("ess_main.Rdata")
ess <- ess  %>%
  mutate(fhighskilled = ifelse(foccupation14=="High-skilled", 1, ifelse(!is.na(foccupation14), 0, NA))) %>%
  mutate(flowskilled = ifelse(foccupation14=="Low-skilled", 1, ifelse(!is.na(foccupation14), 0, NA))) %>%
  mutate(fmodskilled = ifelse(foccupation14=="Medium-skilled", 1, ifelse(!is.na(foccupation14), 0, NA))) %>%
  mutate(mhighskilled = ifelse(moccupation14=="High-skilled", 1, ifelse(!is.na(moccupation14), 0, NA))) %>%
  mutate(mlowskilled = ifelse(moccupation14=="Low-skilled", 1, ifelse(!is.na(moccupation14), 0, NA))) %>%
  mutate(mmodskilled = ifelse(moccupation14=="Medium-skilled", 1, ifelse(!is.na(moccupation14), 0, NA))) %>%
  as.data.frame()

ess <- select(ess, immigration, shock4_1825, shock4_1825comb, fhighskilled, fmodskilled, flowskilled, mhighskilled, mmodskilled, mlowskilled, female, minority, pcitizen, peduc, age, cohort, secondary, eduq_1017, aveedu_1017, educationy, income, lrideology, datalim)

ess$age[ess$age<26] <- NA
ess <- ess[!is.na(ess$age),]
ess$datalim <- NULL

nonNAs <- function(x) {
  as.vector(apply(x, 2, function(x) mean(x, na.rm=TRUE)))}
mean <- nonNAs(ess)
nonNAs <- function(x) {
  as.vector(apply(x, 2, function(x) sd(x, na.rm=TRUE)))}
sd <- nonNAs(ess)
nonNAs <- function(x) {
  as.vector(apply(x, 2, function(x) min(x, na.rm=TRUE)))}
min <- nonNAs(ess)
nonNAs <- function(x) {
  as.vector(apply(x, 2, function(x) max(x, na.rm=TRUE)))}
max <- nonNAs(ess)
nonNAs <- function(x) {
  as.vector(apply(x, 2, function(x) length(which(!is.na(x)))))}
length <- nonNAs(ess)
names <- c("Anti-Immigration Attitude", "Economic Shock (18-25)", "Economic Shock Count (18-25)", "Father High-Skilled", "Father Mod-Skilled", "Father Low-Skilled", "Mother High-Skilled", "Mother Mod-Skilled", "Mother Low-Skilled", "Female", "Minority", "Parental Citzenship", "Parental Education", "Age", "Cohort", "Education 12 years or more", "Educational Equality (10-17)", "Education Level (10-17)", "Education Years", "Income", "Ideology")
table <- data.frame(names, mean, sd, min, max, length)
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA2b.tex", sep = ""))

###################################
##############TABLE A4#############
###################################
##Results when excluding WWII
load("ess_main.Rdata")

ess <- ess[ess$cohort>1935,]
treatment <- c("shock4_1825", "log(shock4_1825comb+1)")
mod <- paste0("immigration ~ ", treatment, baseline, sep="")
label <- c("Economic Shock", "Economic Shock Count")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:2){
  model <- lm(mod[i], data=subset(ess, age > 25 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- label[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  print(counter)
  counter <- counter + 1
}}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA4a.tex", sep = ""))

##Results when those socialized in the 1960s/70s

load("ess_main.Rdata")

ess <- ess[ess$cohort<1935 | ess$cohort>1955,]
mod <- paste0("immigration ~ ", treatment, baseline, sep="")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:2){
  model <- lm(mod[i], data=subset(ess, age > 25 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- label[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  print(counter)
  counter <- counter + 1
}}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA4b.tex", sep = ""))

##Different Definitions of Economic Shock Pre/Post WWII
load("ess_main.Rdata")

##Different definitions for pre/post WWII
treatmentalt <- c("shocka_1825", "log(shocka_1825comb+1)")
mod <- paste0("immigration ~ ", treatmentalt, baseline, sep="")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:2){
  model <- lm(mod[i], data=subset(ess, age > 25 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- label[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  print(counter)
  counter <- counter + 1
}}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA4c.tex", sep = ""))

##4 percent contraction peak to trough

treatmentalt <- c("pt4_1825", "log(pt4_1825comb+1)")
mod <- paste0("immigration ~ ", treatmentalt, baseline, sep="")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:2){
  model <- lm(mod[i], data=subset(ess, age > 25 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- label[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  print(counter)
  counter <- counter + 1
}}
table <- table %>% mutate(lci = V2 - 1.96*V3) %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA4d.tex", sep = ""))


###################################
##############Table A5#############
###################################
##Different definitions of impressionable years.

treatmentalt <- c("shock4_1025", "shock4_1425", "shock4_1625", "shock4_1725", "shock4_1926", "shock4_1821", "shock4_2225", "shock4_2629", "shock4_3033")
treatmentalt <- c(treatmentalt, paste0("log(", treatmentalt, "comb+1)", sep=""))
mod <- paste0("immigration ~ ", treatmentalt, baseline, sep="")
label <- c(rep("Economic Shock", times=9), rep("Economic Shock Count", times=9))
ager <- rep(c("10-25", "14-25", "16-25", "17-25", "19-26", "18-21", "22-25", "26-29", "30-33"), times=2)
ager <- paste0("(", ager, ")", sep="")
label <- paste(label, ager, sep=" ")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:length(mod)){
  model <- lm(mod[i], data=subset(ess, age > 25 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- label[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  print(counter)
  counter <- counter + 1
}}
table <- table %>% mutate(lci = V2 - 1.96*V3) %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table$lci <- gsub(" ", "", table$lci)
table <- table %>% mutate(V3 = paste("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA5.tex", sep = ""))

####################################
##############Table A6##############
####################################
##Using first factor from factor analysis

ess <- select(ess, imm13a, imm13b, imm4, imm5, imm6, shock4_1825, shock4_1825comb, foccupation14, moccupation14, secondary, female, minority, pcitizen, peduc, age, cohort, country, pspwght, essround, aveedu_1017, eduq_1017, c, datalim)
ess <- ess[ess$age>25,]
ess <- na.omit(ess)

fa2 <- factanal(~imm13a+imm13b+imm4+imm5+imm6, data=ess, factors=1, scores="regression")
ess$scores <- fa2$scores

treatment <- c("shock4_1825", "log(shock4_1825comb+1)")
mod <- paste0("scores ~ ", treatment, baseline, sep="")
label <- c("Economic Shock", "Economic Shock Count")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:2){
  model <- lm(mod[i], data=subset(ess, age > 25 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- label[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  print(counter)
  counter <- counter + 1
}}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA6.tex", sep = ""))

####################################
##############Table A7##############
####################################
##Excluding each item separately

load("ess_main.Rdata")

ess$immigration1 <- ess$imm13b + ess$imm4 + ess$imm5 + ess$imm6
ess$immigration2 <- ess$imm13a + ess$imm4 + ess$imm5 + ess$imm6
ess$immigration3 <- ess$imm13a + ess$imm13b + ess$imm5 + ess$imm6
ess$immigration4 <- ess$imm13a + ess$imm13b + ess$imm4 + ess$imm6
ess$immigration5 <- ess$imm13a + ess$imm13b + ess$imm4 + ess$imm5

treatment <- c("shock4_1825", "log(shock4_1825comb+1)")
label <- c("Economic Shock", "Economic Shock Count")
dv <- paste0("immigration", c(1:5), sep="")
table <- data.frame()
counter <- 1

for(v in 1:length(dv)){for(d in 1:2){for(i in 1:2){
  mod <- paste0(dv[v], " ~ ", treatment[i], baseline, sep="")
  model <- lm(mod, data=subset(ess, age > 25 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- label[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  table[counter, 6] <- v
  print(counter)
  counter <- counter + 1
}}}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,7:8)] <- format(round(table[, c(2:3,7:8)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:6)]

table <- table[, -c(6)]

table1 <- table[c(1:4),]
table2 <- table[c(5:8),]
table3 <- table[c(9:12),]
table4 <- table[c(13:16),]
table5 <- table[c(17:20),]

table1 <- xtable(table1, type = "latex", latex.environments = "center", caption = "")
table2 <- xtable(table2, type = "latex", latex.environments = "center", caption = "")
table3 <- xtable(table3, type = "latex", latex.environments = "center", caption = "")
table4 <- xtable(table4, type = "latex", latex.environments = "center", caption = "")
table5 <- xtable(table5, type = "latex", latex.environments = "center", caption = "")

print(table1, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA71.tex", sep = ""))
print(table2, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA72.tex", sep = ""))
print(table3, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA73.tex", sep = ""))
print(table4, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA74.tex", sep = ""))
print(table5, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA75.tex", sep = ""))

###################################
#############Table A8#############
###################################
##Testing items separately
load("ess_main.Rdata")

ess$immigration1 <- ess$imsmetn
ess$immigration2 <- ess$imdfetn
ess$immigration3 <- ess$impcntr
ess$immigration4 <- 10 - ess$imbgeco
ess$immigration5 <- 10 - ess$imueclt
ess$immigration6 <- 10 - ess$imwbcnt

treatment <- c("shock4_1825", "log(shock4_1825comb+1)")
label <- c("Economic Shock", "Economic Shock Count")
dv <- paste0("immigration", c(1:6), sep="")
table <- data.frame()
counter <- 1

for(v in 1:length(dv)){for(d in 1:2){for(i in 1:2){
  mod <- paste0(dv[v], " ~ ", treatment[i], baseline, sep="")
  model <- lm(mod, data=subset(ess, age > 25 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- label[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  table[counter, 6] <- v
  print(counter)
  counter <- counter + 1
}}}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,7:8)] <- format(round(table[, c(2:3,7:8)], 3), nsmall = 3)
table$lci <- gsub(" ", "", table$lci)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:6)]

table$V6 <- NULL

table1 <- table[c(1:4),]
table2 <- table[c(5:8),]
table3 <- table[c(9:12),]
table4 <- table[c(13:16),]
table5 <- table[c(17:20),]
table6 <- table[c(21:24),]
table1 <- xtable(table1, type = "latex", latex.environments = "center", caption = "")
table2 <- xtable(table2, type = "latex", latex.environments = "center", caption = "")
table3 <- xtable(table3, type = "latex", latex.environments = "center", caption = "")
table4 <- xtable(table4, type = "latex", latex.environments = "center", caption = "")
table5 <- xtable(table5, type = "latex", latex.environments = "center", caption = "")
table6 <- xtable(table6, type = "latex", latex.environments = "center", caption = "")
print(table1, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA81.tex", sep = ""))
print(table2, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA82.tex", sep = ""))
print(table3, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA83.tex", sep = ""))
print(table4, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA84.tex", sep = ""))
print(table5, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA85.tex", sep = ""))
print(table6, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA86.tex", sep = ""))

###################################
#############Table A9##############
###################################
##Testing alternative latent measures of immigration attitudes

ess$imm13c <- ess$imm1 + ess$imm2 + ess$imm3
ess$imm13c[ess$imm13c>=0 & ess$imm13c<=2] <- 0
ess$imm13c[ess$imm13c==3] <- 1
ess$imm13d <- ess$imm11 + ess$imm22 + ess$imm33
ess$imm13d[ess$imm13d>=0 & ess$imm13d<=2] <- 0
ess$imm13d[ess$imm13d==3] <- 1
ess$immigration7 <- ess$imm13c + ess$imm13d + ess$imm4 + ess$imm5 + ess$imm6
ess$immigration8 <- ess$imm11 + ess$imm22 + ess$imm33 + ess$imm4 + ess$imm5 + ess$imm6
ess$imm111 <- scale(ess$immigration1)[, 1]
ess$imm222 <- scale(ess$immigration2)[, 1]
ess$imm333 <- scale(ess$immigration3)[, 1]
ess$imm444 <- scale(ess$immigration4)[, 1]
ess$imm555 <- scale(ess$immigration5)[, 1]
ess$imm666 <- scale(ess$immigration6)[, 1]
ess$immigration9 <- (ess$imm111 + ess$imm222 + ess$imm333 + ess$imm444 + ess$imm555 + ess$imm666)/6

treatment <- c("shock4_1825", "log(shock4_1825comb+1)")
label <- c("Economic Shock", "Economic Shock Count")
dv <- paste0("immigration", c(7:9), sep="")
table <- data.frame()
counter <- 1

for(v in 1:length(dv)){for(d in 1:2){for(i in 1:2){
  mod <- paste0(dv[v], " ~ ", treatment[i], baseline, sep="")
  model <- lm(mod, data=subset(ess, age > 25 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- label[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  table[counter, 6] <- v
  print(counter)
  counter <- counter + 1
}}}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,7:8)] <- format(round(table[, c(2:3,7:8)], 3), nsmall = 3)
table$lci <- gsub(" ", "", table$lci)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]

table1 <- table[c(1:4),]
table2 <- table[c(5:8),]
table3 <- table[c(9:12),]
table1 <- xtable(table1, type = "latex", latex.environments = "center", caption = "")
table2 <- xtable(table2, type = "latex", latex.environments = "center", caption = "")
table3 <- xtable(table3, type = "latex", latex.environments = "center", caption = "")
print(table1, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA91.tex", sep = ""))
print(table2, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA92.tex", sep = ""))
print(table3, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA93.tex", sep = ""))

  
###################################
#############Table A10#############
###################################
##Country-cohort analysis
load("ess_main.Rdata")
completeFun <- function(data, desiredCols) {
  completeVec <- complete.cases(data[, desiredCols])
  return(data[completeVec, ])
}

library("collapse")

ess <- ess[ess$age>25,]
ess <- select(ess, country, pspwght, immigration, cohort, c, shock4_1825, shock4_1825comb, datalim, aveedu_1017, eduq_1017)
ess <- completeFun(ess, c("immigration"))
ess$id <- 1
ess$id1 <- ave(ess$id, ess$c, FUN = sum)
ess <- collap(ess, ~ country + cohort, custom = list(fsum = c("id"), fmean = c("immigration", "datalim", "shock4_1825", "shock4_1825comb", "id1", "aveedu_1017", "eduq_1017")), w = ~ pspwght)

mod <- c("shock4_1825", "log(shock4_1825comb+1)")
mod1 <- paste0("immigration ~ ", mod, " + aveedu_1017 + eduq_1017 + factor(country)", sep="")
mod2 <- paste0("immigration ~ ", mod, " + aveedu_1017 + eduq_1017 + factor(cohort) + factor(country)", sep="")
name <- c("Economic Shock", "Economic Shock Count")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:2){
  model <- lm(mod1[i], data=subset(ess, datalim<d), weights=id)
  model <- robustify(model)
  table[counter, 1] <- name[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- "No"
  table[counter, 6] <- ifelse(d==1, "Restricted", "Full")

  counter <- counter + 1
  model <- lm(mod2[i], data=subset(ess, datalim<d), weights=id)
  model <- robustify(model)
  table[counter, 1] <- name[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- "Yes"
  table[counter, 6] <- ifelse(d==1, "Restricted", "Full")
  print(counter)
  counter <- counter + 1
}}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,7:8)] <- format(round(table[, c(2:3,7:8)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:6)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA10.tex", sep = ""))

###################################
#############Table A11#############
###################################
##Excluding each country
load("ess_main.Rdata")

cc1 <- levels(as.factor(ess$country[ess$datalim==0]))
cc2 <- gsub("(?<=\\b)([a-z])", "\\U\\1", tolower(cc1), perl=TRUE)
mod <- paste0("immigration ~ shock4_1825", baseline, sep="")
table <- data.frame()
counter <- 1

for(i in 1:length(cc1)){
  ccc1 <- cc1[i]
  ccc2 <- cc2[i]
  model <- lm(mod, data=subset(ess, age > 25 & datalim<1 & country!=ccc1), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- paste0("Excluding ", ccc2, sep="")
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  print(counter)
  counter <- counter + 1
}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,5:6)] <- format(round(table[, c(2:3,5:6)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:4)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA11a.tex", sep = ""))

cc1 <- levels(as.factor(ess$country))
cc2 <- gsub("(?<=\\b)([a-z])", "\\U\\1", tolower(cc1), perl=TRUE)
mod <- paste0("immigration ~ shock4_1825", baseline, sep="")
table <- data.frame()
counter <- 1

for(i in 1:length(cc1)){
  ccc1 <- cc1[i]
  ccc2 <- cc2[i]
  model <- lm(mod, data=subset(ess, age > 25 & datalim<2 & country!=ccc1), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- paste0("Excluding ", ccc2, sep="")
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  print(counter)
  counter <- counter + 1
}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,5:6)] <- format(round(table[, c(2:3,5:6)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:4)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA11b.tex", sep = ""))

##################################
#############Table A12############
##################################
##Excluding each round

r <- c(1:8)
mod <- paste0("immigration ~ shock4_1825", baseline, sep="")
table <- data.frame()
counter <- 1

for(i in 1:length(r)){
  rr <- r[i]
  model <- lm(mod, data=subset(ess, age > 25 & datalim<1 & essround!=rr), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- paste0("Excluding Round ", rr, sep="")
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  print(counter)
  counter <- counter + 1
}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,5:6)] <- format(round(table[, c(2:3,5:6)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:4)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA12a.tex", sep = ""))

table <- data.frame()
counter <- 1

for(i in 1:length(r)){
  rr <- r[i]
  model <- lm(mod, data=subset(ess, age > 25 & datalim<2 & essround!=rr), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- paste0("Excluding Round ", rr, sep="")
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  print(counter)
  counter <- counter + 1
}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,5:6)] <- format(round(table[, c(2:3,5:6)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:4)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA12b.tex", sep = ""))

#####################################
##############TABLE A13##############
#####################################

baselinealt <- " + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + dem_1017 + eqdr_1017 + gdp_1017 + factor(cohort) + age + I(age^2) + factor(country)*age + factor(essround)"

treatment <- c("shock4_1825", "log(shock4_1825comb+1)")
mod <- paste0("immigration ~ ", treatment, baselinealt, sep="")
label <- c("Economic Shock", "Economic Shock Count")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:2){
	model <- lm(mod[i], data=subset(ess, age > 25 & datalim<d), weights=pspwght)
	model <- robustify(model, cluster=c)
	table[counter, 1] <- label[i]
	table[counter, 2] <- summary(model)$coefficients[2,1]
	table[counter, 3] <- summary(model)$coefficients[2,2]
	table[counter, 4] <- nobs(model)
	table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
	print(counter)
	counter <- counter + 1
	}}
table <- table %>% mutate(lci = V2 - 1.96*V3) %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA13.tex", sep = ""))


##################################
#############Table A15############
##################################

baseline <- " + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround)"

mod <- c("shock4_1825", "log(shock4_1825comb+1)")
mod <- paste0("immigration ~ ", mod, baseline, sep="")
name <- c("Economic Shock", "Economic Shock Count")
table <- data.frame()
counter <- 1

for(d in 1:2){for(i in 1:2){
  model <- lm(mod[i], data=subset(ess, age > 17 & datalim<d), weights=pspwght)
  model <- robustify(model, cluster=c)
  table[counter, 1] <- name[i]
  table[counter, 2] <- summary(model)$coefficients[2,1]
  table[counter, 3] <- summary(model)$coefficients[2,2]
  table[counter, 4] <- nobs(model)
  table[counter, 5] <- ifelse(d==1, "Restricted", "Full")
  print(counter)
  counter <- counter + 1
}}
table <- table %>% mutate(lci = V2 - 1.96*V3)
table <- table %>% mutate(uci = V2 + 1.96*V3)
table[, c(2:3,6:7)] <- format(round(table[, c(2:3,6:7)], 3), nsmall = 3)
table <- table %>% mutate(V3 = paste0("(", lci, ", ", uci, ")", sep=""))
table <- tab <- table[, c(1:5)]
table <- xtable(table, type = "latex", latex.environments = "center", caption = "")
print(table, include.rownames=FALSE, include.colnames=FALSE,  sanitize.text.function=identity, only.contents=TRUE,  hline.after = NULL, file=paste(location, "paper/tables/tableA15.tex", sep = ""))

#####################################
##############Figure A2##############
#####################################
load("ess_main.Rdata")

baseline <- " + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround)"

mod <- paste0("immigration ~ ", "shock4_1825", baseline, sep="")
model <- lm(mod, data=subset(ess, age > 25 & datalim<1), weights=pspwght)
model <- robustify(model, cluster=c)

v <- c("Economic Shock (18-25)", "12 years or more of edu.", "Father: Low-skilled", "Father: Mod-skilled", "Mother: Low-skilled", "Mother: Mod-skilled", "Female", "Minority", "Parental Citizenship", "Parental Education")
table <- data.frame()
counter <- 1
for(i in 2:11){
  table[counter, 1] <- v[counter]
  table[counter, 2] <- summary(model)$coefficients[i,1]
  table[counter, 3] <- summary(model)$coefficients[i,2]
  counter <- counter + 1
}
V1 <- c("Father: High-skilled", "Mother: High-skilled")
V2 <- c(0, 0)
V3 <- c(0, 0)
extra <- data.frame(V1, V2, V3)
table <- rbind(table, extra)
table <- table[c(1,3:4,11,5:6,12,7:10,2),]
table[12,2] <- table[12,2]/6
table[12,3] <- table[12,3]/6
table <- table %>% mutate(lci95 = V2 - 1.96*V3) %>% mutate(uci95 = V2 + 1.96*V3) %>% mutate(lci90 = V2 - 1.645*V3) %>% mutate(uci90 = V2 + 1.645*V3)

setwd(paste(location, "paper/figures", sep = ""))
pdf("figureA2.pdf", height=6, width=12)
par(mgp=c(3,1,0))
par(mar=c(5,5,0,.25))
counter <- 11
plot(x = NA, xlim = c(-0.30, 0.25), ylim = c(-0.25,12), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=2, at=c(11:0), labels=table$V1, pos=-0.24, lty=1, las=2)
axis(side=1, at=round(seq(from=-0.20, to=0.30, by=0.10),3), lty=1, las=0)
abline(v = 0, col = "grey", lty = 2)
mtext(side=1, at=0, text="Estimated Coefficient", line=3, cex=1.25)
for(i in 1:nrow(table)){
  points(x = table[i, 2], y = counter, pch = 19, cex=1)
  segments(x0 = table[i, 4], x1 = table[i, 5], y0 = counter, cex=1.25)
  segments(x0 = table[i, 6], x1 = table[i, 7], y0 = counter, lwd=2)
  counter <- counter - 1	
}
dev.off()

#####################################
##############Figure A3##############
#####################################

table <- as.data.frame(matrix(NA, nrow=16, ncol=8))
colnames(table) <- c("age", "treatment", "sample", "est", "lci1", "uci1", "lci5", "uci5")
counter <- 1
a <- seq(from=30, to=65, by=5)
baseline <- "*age + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + I(age^2) + factor(country) + factor(essround)"
treatment <- c("shock4_1825")
mod <- paste0("immigration ~ ", treatment, baseline, sep="")


for(d in 1:2){
	model <- lm(mod, data=subset(ess, age > 25 & datalim<d), weights=pspwght)
	model <- robustify(model, cluster=c)
	mod1 <- slopes(model, variables=treatment, newdata = datagrid(age=c(seq(from=30, to=65, by=5))))
	for(p in 1:nrow(mod1)){
		table[counter, 1] <- a[p]
		table[counter, 2] <- treatment
		table[counter, 3] <- ifelse(d==1, "Restricted", "Full")	
		table[counter, 4] <- mod1[p,4]
		table[counter, 5] <- mod1[p,4] - 1.645*mod1[p,5]
		table[counter, 6] <- mod1[p,4] + 1.645*mod1[p,5]
		table[counter, 7] <- mod1[p,4] - 1.96*mod1[p,5]
		table[counter, 8] <- mod1[p,4] + 1.96*mod1[p,5]
		print(counter)
		counter <- counter + 1
	}
}


setwd(paste(location, "paper/figures", sep = ""))
pdf("figureA3.pdf", height=5, width=10)
par(mar=c(5,5,0,.25))
plot(x = NA, ylim = c(-0.10, 0.27), xlim = c(29,65), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=1, at=seq(from=30,to=65,by=5), pos=-0.105, lty=1, las=1, cex.axis=1.15)
axis(side=2, at=c(-0.10, -0.05, 0.00, 0.05, 0.10, 0.15, 0.20, 0.25), lty=1, las=1, cex.axis=1.15)
abline(h = 0, col = "grey", lty = 2)
mtext(side=1, at=47.5, text="Age", line=3.5, cex=1.25)
mtext(side=2, text="Marginal Effect of Economic Shock", line=3.5, cex=1.25)
counter <- 30
j <- 0.5
for(i in 1:8){
  points(y = table[i,4], x = counter+j, pch = 19, cex=1.25,)
  segments(y0 = table[i,7], y1 = table[i,8], x0 = counter+j, cex=1.25)
  segments(y0 = table[i,5], y1 = table[i,6], x0 = counter+j, lwd=2)
  points(y = table[i+8,4], x = counter-j, pch = 19, cex=1.25, col="gray")
  segments(y0 = table[i+8,7], y1 = table[i+8,8], x0 = counter-j, cex=1.25, col="gray")
  segments(y0 = table[i+8,5], y1 = table[i+8,6], x0 = counter-j, lwd=2, col="gray")
  counter <- counter + 5
}
dev.off()


#####################################
##############Figure A4##############
#####################################

baseline <- " + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround)"

mod <- paste0("immigration ~ shock4_1825 + shockoutside", baseline, sep="")
label <- c("Economic Shock 18-25", "Economic Shock Outside Impressionable Years")
table <- data.frame()
counter <- 1

for(d in 1:2){
	model <- lm(mod, data=subset(ess, age > 25 & datalim<d), weights=pspwght)
	model <- robustify(model, cluster=c)
	for(i in 1:2){
		table[counter, 1] <- summary(model)$coefficients[i+1,1]
		table[counter, 2] <- summary(model)$coefficients[i+1,1] - 1.96*summary(model)$coefficients[i+1,2]
		table[counter, 3] <- summary(model)$coefficients[i+1,1] + 1.96*summary(model)$coefficients[i+1,2]
		table[counter, 4] <- summary(model)$coefficients[i+1,1] - 1.645*summary(model)$coefficients[i+1,2]
		table[counter, 5] <- summary(model)$coefficients[i+1,1] + 1.645*summary(model)$coefficients[i+1,2]
		table[counter, 6] <- nobs(model)
		table[counter, 7] <- ifelse(d==1, "Restricted", "Full")
		table[counter, 8] <- label[i]
		print(counter)
		counter <- counter + 1
	}}
	
table1 <- table[c(1,2),]
table2 <- table[c(3,4),]

setwd(paste(location, "paper/figures", sep = ""))
pdf("figureA4.pdf", height=4, width=9) 
par(mfrow=c(1,1))
par(mgp=c(3,1,0))
par(mar=c(5,1,1,0))
plot(x = NA, xlim = c(-0.05, 0.205), ylim = c(-.5,11.5), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
mtext("Estimated Coefficient", side=1, line=2.5, cex=1.05)
axis(side=1, at=c(-0.05, 0.00, 0.05, 0.10, 0.15, 0.20), lty=1, las=0)
abline(v = 0, col = "grey", lty = 2)
counter <- 1
for(i in 1:length(table1)){
	points(x = table2[i,1], y = counter, pch = 19, cex=0.80, col="gray")
	segments(x0 = as.numeric(table2[i,2]), x1 = table2[i,3], y0 = counter, col="gray")
	segments(x0 = as.numeric(table2[i,4]), x1 = table2[i,5], y0 = counter, lwd=2, col="gray")
	text(x = table2[i,1], y = counter+1, table2[i,8], cex=0.80, col="gray")

	points(x = table1[i,1], y = counter+6, pch = 19, cex=0.80)
	segments(x0 = as.numeric(table1[i,2]), x1 = table1[i,3], y0 = counter+6, cex=0.80)
	segments(x0 = as.numeric(table1[i,4]), x1 = table1[i,5], y0 = counter+6, lwd=2)
	text(x = table1[i,1], y = counter+7, table1[i,8], cex=0.80)
	counter <- counter + 1.5
}
dev.off()


###################################
#############Figure A5#############
###################################

t <- c(1:6)
t <- paste0("log(shock", t, "_1825comb+1)", sep="")
counter <- 1
table <- data.frame()

for(d in 1:2){for(i in 1:length(t)){
	mod <- paste0("immigration ~ ", t[i], baseline, sep="")
	model <- lm(mod, data=subset(ess, age > 25 & datalim<d), weights=pspwght)
	model <- robustify(model, cluster=c)
	table[counter, 1] <- summary(model)$coefficients[2,1]
	table[counter, 2] <- summary(model)$coefficients[2,1] - 1.96*summary(model)$coefficients[2,2]
	table[counter, 3] <- summary(model)$coefficients[2,1] + 1.96*summary(model)$coefficients[2,2]
	table[counter, 4] <- summary(model)$coefficients[2,1] - 1.645*summary(model)$coefficients[2,2]
	table[counter, 5] <- summary(model)$coefficients[2,1] + 1.645*summary(model)$coefficients[2,2]
	table[counter, 6] <- t[i]
	table[counter, 7] <- ifelse(d==1, "Restricted", "Full")
	print(counter)
	counter <- counter + 1
	}}

tab1 <- table[1:6,]
tab2 <- table[7:12,]

setwd(paste(location, "paper/figures", sep = ""))
pdf("figureA5.pdf", height=5, width=10)
par(mfrow=c(1,1))
par(mgp=c(3,1,0))
par(mar=c(5,6,0,.25))
plot(x = NA, xlim = c(-0.05, 0.20), ylim = c(-.25,5), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=2, at=c(0,1,2,3,4,5), labels=c("1 Percent","2 Percent","3 Percent","4 Percent","5 Percent","6 Percent"), pos=-0.055, lty=1, las=1)
axis(side=1, at=c(-0.05, 0.00, 0.05, 0.10, 0.15, 0.20), lty=1, las=1)
abline(v = 0, col = "grey", lty = 2)
mtext(side=2, text="Contraction Threshold", line=5, cex=1.25)
mtext(side=1, text="Coefficient on Economic Shock Count", line=3, cex=1.25)
counter <- 0
for(i in 1:nrow(table)){
	points(x = tab1[i,1], y = counter+0.10, pch = 19, cex=1.25)
	segments(x0 = tab1[i,2], x1 = tab1[i,3], y0 = counter+0.10, cex=1.25)
	segments(x0 = tab1[i,4], x1 = tab1[i,5], y0 = counter+0.10, lwd=2)
	points(x = tab2[i,1], y = counter-0.10, pch = 19, cex=1.25, col="gray")
	segments(x0 = tab2[i,2], x1 = tab2[i,3], y0 = counter-0.10, cex=1.25, col="gray")
	segments(x0 = tab2[i,4], x1 = tab2[i,5], y0 = counter-0.10, lwd=2, col="gray")
	counter <- counter + 1
}
dev.off()


#####################################
##############Figure A6##############
#####################################

table <- as.data.frame(matrix(NA, nrow=2, ncol=6))
colnames(table) <- c("sample", "est", "lci1", "uci1", "lci5", "uci5")
baseline <- " + secondary + factor(foccupation14) + factor(moccupation14) + female + minority + pcitizen + peduc + aveedu_1017 + eduq_1017 + factor(cohort) + age + I(age^2) + factor(country) + factor(essround)"
treatment <- c("shock4_1825")
mod <- paste0("racial ~ ", treatment, baseline, sep="")

for(d in 1:2){
	model <- glm(mod, data=subset(ess, age > 25 & datalim<d), family=binomial(link='logit'), weights=pspwght)
	model <- robustify(model, cluster=c)
	mod1 <- avg_slopes(model, variables="shock4_1825")
	for(p in 1:nrow(mod1)){
		table[d, 1] <- ifelse(d==1, "Restricted", "Full")	
		table[d, 2] <- mod1[p,3]
		table[d, 3] <- mod1[p,3] - 1.645*mod1[p,4]
		table[d, 4] <- mod1[p,3] + 1.645*mod1[p,4]
		table[d, 5] <- mod1[p,3] - 1.96*mod1[p,4]
		table[d, 6] <- mod1[p,3] + 1.96*mod1[p,4]
	}
}


setwd(paste(location, "paper/figures", sep = ""))
pdf("figureA6.pdf", height=5, width=10)
par(mfrow=c(1,1))
par(mgp=c(3,1,0))
par(mar=c(5,1,0,.25))
plot(x = NA, xlim = c(-0.01, 0.04), ylim = c(.75,2.25), bty = "n", yaxt = "n", xaxt="n", ylab = "", xlab = "")
axis(side=1, at=c(-0.01, 0.00, 0.01, 0.02, 0.03, 0.04), lty=1, las=1)
abline(v = 0, col = "grey", lty = 2)
mtext(side=1, text="Marginal Effect of Economic Shock", line=3, cex=1.25)
j <- 0.05
for(i in nrow(table):1){
  points(x = table[i,2], y = i, pch = 19, cex=1.25)
  segments(x0 = table[i,5], x1 = table[i,6], y0 = i, cex=1.25)
  segments(x0 = table[i,3], x1 = table[i,4], y0 = i, lwd=2)
  text(x = table[i,2], y = i+j, paste0("Economic Shock (18-25) ", table[i,1], sep=""), cex=1.25)
}
dev.off()







